
###############################################################
# Figures 3 and 4: Allocation of political attention, 2001-2014
###############################################################

# Programme:  cfsp-analysis03a-wpRel.r
# Author:     Frank Haege, Department of Politics and Administration, University of Limerick
# Contact:		frank.haege@ul.ie

# Description
#############
# This script creates Figure 3 (attention to horizontal issues working parties) and 
# Figure 4 (attention to geographical areas working parties) of the article. 


# Remove all objects
rm(list = ls(all = TRUE))

# Load libraries
library(data.table)
library(ggplot2)
library(ggpubr)
library(mgcv)

# Set working directory
setwd("C:\\Users\\Frank\\Dropbox\\European Union\\CFSP attention\\EUP submission")

# Set options
options(stringsAsFactors=F)

# Load data
load("Data management\\cfsp-management02-semesterData.RData")

# Drop general meeting variables
df = df[!(df$id=="C00" | df$id=="D00" | df$id=="C51"), ]
# These categories include meetings that could not be unambiguously allocated to any particular working party

# Remove coordinating working parties
df = df[df$coord!=1, ]
# Coordinating working parties do not belong to any particular sub-area

# Restrict data to CFSP (without CSDP)
df = df[df$area=="Horizontal" | df$area=="Regional", ]

# Remove all working party variables
df = subset(df, select=-c(wppolicy, nowp, coord, comp, diff.a, dur.lag, diff.r))
df = df[order(df$yearsem, df$area, df$id), ]

# Calculate share of meetings by semester
df = data.table(df)
df = df[, perc:=prop.table(dur)*100, by="yearsem"]
df = data.frame(df)


# Horizontal issues working parties
###################################

# Loop through all horizontal working parties
wpname.list = names(table(df$wpname[df$area=="Horizontal"]))
plot.list = list()
for (i in 1:11) {
  j = wpname.list[[i]]
  dta = df[df$wpname==j, ]
  
  # Estimate GAM of number of meetings with interaction of time by Lisbon Treaty
  m = gam(data=dta, formula = perc ~ s(time) + s(time, by=lisbonsem) + lisbonsem)
  #print(summary(m))
  
  # Predict fitted values and standard errors
  pv = data.frame(predict(m, se.fit=T, type="response"))
  
  # Create dataframe for fitted values and confidence intervals
  fitted = data.frame(perc=dta$perc, yearsem=dta$yearsem, time=dta$time, lisbonsem=dta$lisbonsem, 
                      pv, wpname=dta$wpname)
  fitted$ci.low = fitted$fit-1.96*fitted$se.fit
  fitted$ci.high = fitted$fit+1.96*fitted$se.fit
  
  # Plot duration of meetings over time and overlay predicted values
  p = ggplot() +
    geom_point(data=fitted, aes(x=time, y=perc)) +
    ggtitle(j) +
    geom_line(data=fitted[fitted$lisbonsem==0, ], aes(x=time, y=fit)) +
    geom_line(data=fitted[fitted$lisbonsem==1, ], aes(x=time, y=fit)) +
    geom_ribbon(data=fitted[fitted$lisbonsem==0, ],aes(x=time, ymin=ci.low,ymax=ci.high),alpha=0.3) +
    geom_ribbon(data=fitted[fitted$lisbonsem==1, ],aes(x=time, ymin=ci.low,ymax=ci.high),alpha=0.3) +
    #ylim(0, 250) + 
    ylab("Share of meetings") +
    xlab("Year") +
    geom_vline(xintercept=18.5, linetype="dashed") +
    scale_x_continuous(breaks=c(3, 11, 19, 27), labels=c("2002", "2006", "2010", "2014")) +
    theme_bw(base_size=10) +
    theme(panel.grid.major = element_blank(), 
          panel.grid.minor = element_blank())
  plot.list[[i]] = p
}
p01 = plot.list[[1]]
p02 = plot.list[[2]]
p03 = plot.list[[3]]
p04 = plot.list[[4]]
p05 = plot.list[[5]]
p06 = plot.list[[6]]
p07 = plot.list[[7]]
p08 = plot.list[[8]]
p09 = plot.list[[9]]
p10 = plot.list[[10]]
p11 = plot.list[[11]]
rm(i, j, m, p, plot.list, wpname.list, fitted, dta, pv)

# Combine and export plots
plots = ggarrange(p01, p02, p03, p04, p05, p06, p07, p08, p09, p10, p11,
                  ncol=2, nrow=6)
plots
ggexport(plots, filename="Data analysis\\Figures\\cfsp-analysis03a-wpRel-horizontal.png", width=500, height=1000)


# Geographical areas working parties
####################################

# Loop through all regional working parties
wpname.list = names(table(df$wpname[df$area=="Regional"]))
plot.list = list()
for (i in 1:9) {
  j = wpname.list[[i]]
  dta = df[df$wpname==j, ]
  
  # Estimate GAM of number of meetings with interaction of time by Lisbon Treaty
  m = gam(data=dta, formula = perc ~ s(time) + s(time, by=lisbonsem) + lisbonsem)
  #print(summary(m))
  
  # Predict fitted values and standard errors
  pv = data.frame(predict(m, se.fit=T, type="response"))
  
  # Create dataframe for fitted values and confidence intervals
  fitted = data.frame(perc=dta$perc, yearsem=dta$yearsem, time=dta$time, lisbonsem=dta$lisbonsem, 
                      pv, wpname=dta$wpname)
  fitted$ci.low = fitted$fit-1.96*fitted$se.fit
  fitted$ci.high = fitted$fit+1.96*fitted$se.fit
  
  # Plot duration of meetings over time and overlay predicted values
  p = ggplot() +
    geom_point(data=fitted, aes(x=time, y=perc)) +
    ggtitle(j) +
    geom_line(data=fitted[fitted$lisbonsem==0, ], aes(x=time, y=fit)) +
    geom_line(data=fitted[fitted$lisbonsem==1, ], aes(x=time, y=fit)) +
    geom_ribbon(data=fitted[fitted$lisbonsem==0, ],aes(x=time, ymin=ci.low,ymax=ci.high),alpha=0.3) +
    geom_ribbon(data=fitted[fitted$lisbonsem==1, ],aes(x=time, ymin=ci.low,ymax=ci.high),alpha=0.3) +
    #ylim(0, 250) + 
    ylab("Share of meetings") +
    xlab("Year") +
    geom_vline(xintercept=18.5, linetype="dashed") +
    scale_x_continuous(breaks=c(3, 11, 19, 27), labels=c("2002", "2006", "2010", "2014")) +
    theme_bw(base_size=10) +
    theme(panel.grid.major = element_blank(), 
          panel.grid.minor = element_blank())
  plot.list[[i]] = p
}
p01 = plot.list[[1]]
p02 = plot.list[[2]]
p03 = plot.list[[3]]
p04 = plot.list[[4]]
p05 = plot.list[[5]]
p06 = plot.list[[6]]
p07 = plot.list[[7]]
p08 = plot.list[[8]]
p09 = plot.list[[9]]
#rm(i, j, m, p, plot.list, wpname.list, fitted, dta, pv)

# Combine and export plots
plots = ggarrange(p01, p02, p03, p04, p05, p06, p07, p08, p09,
                  ncol=2, nrow=5)
plots
ggexport(plots, filename="Data analysis\\Figures\\cfsp-analysis03a-wpRel-regional.png", width=500, height=865)


